      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
      double precision E, v, lambda, mu
      double precision misesEqualStress, yield0, yieldn, yieldn1, hydroStaticPressure, H
      double precision equalPlasticStrain,dEqualPlasticStrain
      INTEGER i,j
      double precision, DIMENSION(NTENS)::flow, deviatoricStress

      INTEGER,save::cn=0
      cn=cn+1

C      user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
C      and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
      E=PROPS(1)
      v=PROPS(2)     
      yield0=PROPS(3)

      lambda=E*v/((1+v)*(1-2*v))
      mu=E/(2*(1+v))
      H=(PROPS(5)-PROPS(3))/(PROPS(6)-PROPS(4))

C     Calculate the material Jacobian matrix
      do i=1,NTENS
         do j=1,NTENS
            DDSDDE(i,j)=0
         enddo
      enddo
      do i=1,NDI
         do j=1,NDI
            DDSDDE(i,j)=lambda
         enddo
         DDSDDE(i,i)=lambda+2*mu
         DDSDDE(i+NDI,i+NDI)=mu
      enddo

C     Calculate stress-try      
      do i=1,NTENS
         do j=1,NTENS
            STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
         enddo
      enddo
      
C     Calculate hydroStatic Pressure
      hydroStaticPressure=0
      do i=1,NDI
         hydroStaticPressure=hydroStaticPressure+STRESS(i)
      enddo
      hydroStaticPressure=hydroStaticPressure/3

C     Calculate the deviatoric stress-try 
      do i=1,NDI
         deviatoricStress(i)=STRESS(i)-hydroStaticPressure
      enddo
      do i=NDI+1,NTENS
         deviatoricStress(i)=STRESS(i)
      enddo
      
C     Calculate  Misess equal stress
      misesEqualStress=0
      do i=1,NDI
         misesEqualStress=misesEqualStress+deviatoricStress(i)**2
      enddo
      do i=1+NDI,NTENS
         misesEqualStress=misesEqualStress+2*deviatoricStress(i)**2
      enddo

      misesEqualStress=sqrt(1.5*misesEqualStress)

C     get  the nth equalPlasticStrain from STATEV
      equalPlasticStrain=STATEV(1)

C     Calculate  the nth yield
      yieldn=yield0+H*(equalPlasticStrain)

C     Calculate  the nth yield

      if (misesEqualStress>yieldn) then     
C     Calculate dEqualPlasticStrain          
         dEqualPlasticStrain=(misesEqualStress-yieldn)/(H+3*mu)

C     Calculate the deviatoric stress-(n+1) 
         do i=1,NTENS
            deviatoricStress(i)=deviatoricStress(i)-3*mu*dEqualPlasticStrain*deviatoricStress(i)/misesEqualStress
         enddo
C     Calculate the  stress-(n+1)          
         do i=1,NDI
            STRESS(i)=deviatoricStress(i)+hydroStaticPressure
         enddo
         do i=NDI+1,NTENS
            STRESS(i)=deviatoricStress(i)
         enddo

      else       
         dEqualPlasticStrain=0
      end if

      do i=1,NTENS
         STRAN(i)=STRAN(i)+DSTRAN(i)
      enddo

     

      STATEV(1)=equalPlasticStrain+dEqualPlasticStrain
      RETURN
      END